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Abstract 

A method for photoacoustic tomography is presented that uses circular integrals of the 
acoustic wave for the reconstruction of a three-dimensional image. Image reconstruction is a 
two-step process: In the first step data from a stack of circular integrating detectors are used 
to reconstruct the circular projection of the source distribution. In the second step the inverse 
circular Radon transform is applied. In this article we establish inversion formulas for the first 
step, which involves an inverse problem for the axially symmetric wave equation. Numerical 
results are presented that show the validity and robustness of the resulting algorithm. 

Keywords. Radon transform; Photoacoustic tomography; photoacoustic microscopy; 
Hankel transform; image reconstruction; integrating detectors; axially symmetric; wave equa- 
tion; 

AMS classifications. 44A12, 65R32, 35L05, 92C55. 

1 Introduction 

The principle of Photoacoustic tomography (PAT), also called Thermoacoustic tomography, is 
based on the excitation of high bandwidth acoustic waves with pulsed non-ionizing electromag- 
netic energy inside tissue [9, 17, 23, 25, 28]. PAT presents a hybrid imaging technique that 
combines the advantages of optical (high contrast) and ultrasound imaging (high resolution). It 
has proven great potential for important medical applications including cancer diagnostics [16, 21] 
and imaging of vasculature [7, 15]. 

The common approach uses small conventional piezoelectric transducers that approximate 
point detectors to measure the acoustic pressure [28] . Reconstruction algorithms which are based 
on the point detector assumption yield images with a spatial resolution that is essentially limited 
by the size of the used piezoelectric transducers [27] . The size of the detector can in principle be 
reduced, but only at the cost of also reducing the signal-to-noise ratio. In order to overcome this 
limitation large size planar or linear integrating detectors have been proposed in [3, 12]. Line 
shaped detectors integrate the acoustic pressure over its length and can be implemented by a 
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Mach-Zehnder [22] or a Fabry-Perot interferometer [10], which naturaUy integrate the acoustic 
pressure over the length of a laser beam. 

A drawback of linear integrating detectors is that because of attenuation parts of the detector 
which are distant from the object may be less influenced by the pressure wave. Moreover, linear 
integrating detectors do not provide a compact experimental buildup. In [32] so called circular 
integrating detectors where introduced, which integrate the acoustic pressure over circles. A 
circular integrating detector can, similar to a linear integrating detector, be implemented by an 
interferometer where the laser beam is guided along a circle in an optical fiber. It is free of any 
aperture effect and can provide a uniformly and high resolution throughout the image area. Since 
it is possible to fabricate optical fibers out of materials which have nearly the same acoustical 
density like the fluid in which they are contained [10] no shadowing effects due to the circular 
integrating detectors are expected. The use of circular integrating detectors has been proposed 
independently in [31]. However, their study was limited to two spatial dimensions, where the 
circular shaped detector can be used as virtual point detector. 

In [32] we derived an inversion formula based on the expansion of the involved functions in 
special basis-functions. This formula, is hard to implement directly due to possible division by 
a zero (see Subsection 2.2). In this article we prove two novel exact scries solutions that allow 
for stable implementation. As a byproduct we obtain a novel inversion formula for PAT using 
point-like detectors on a cylindrical recording surface (see Remark 3.6). 

The outline of this article is as follows. In Section 2 we review PAT with circular integrating 
detectors and recall the main results of [32] . In the following Sections we will be concerned with 
derivation and implementation of the novel stable inversion formulas. 
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2 Photoacoustic Tomography with Circular Integrating De- 
tectors 

In PAT an acoustic pressure wave p(x, t), inside an object, is generated by a pulse of non-ionizing 
electromagnetic radiation. In the case of spatially constant sound speed, the induced acoustic 
pressure satisfies the initial value problem [25, 28] 

dM^,t) = Ap{^,t), (x,t) €]R3 X (0,(X3), (2.1) 

p(x,0) = /(x), xeM-\ (2.2) 

atp(x, 0) = , X e M-\ (2.3) 

PAT is concerned with recovering the initial pressure from measurements of p outside the support 
off. 

In [32] it is proposed to measure the acoustic signals with a stack of parallel circles which 
is rotated around a single axis, see Figure 1. In such a situation, three-dimensional imaging 
involves the inversion of the classical circular Radon transform and the inversion of a reduced 
wave equation, as outlined in the following. 

Throughout this article it is assumed that / is smooth and supported in the cylinder -8^(0) x K, 
where i? is a fixed positive number. Let p{x,t) denote the unique solution of (2.1)-(2.3) and, for 
a & S^, define 



1 /■^'^ 

Pa{z,r,t) := — p{^„{z,r,a),t)da, (z, r, t) € R x (0, ( 

271- Jo 



oo)2, (2.4) 



where 



1 f^^ 

— J^ f{^„{z,r,a))da, (z, r) e M x (0, oo) , (2.5) 



^a{z, r, a) = Ra + {r cos{a),r sin{a), z) , (z, r, a) e K x (0, oo) x [0, 27r] . 

The stack of circular integrating detectors measures 

G'^(z,t) := P^(z,rdot,i), (cr, z, i) e 5^ x M x (0, oo) , (2.6) 

with rdet > denoting the fixed radius of the detectors. 

The goal is to recover the unknown initial data / from measured data {Ga)crGS^- 

2.1 Two Stage Reconstruction 

Reconstruction with circular integrating detectors is based on the following reduction to the axial 
symmetric wave equation: 

Proposition 2.1. Let f € C§°{Br{0) x M) and define P„ and F„, a e by (24), (2.5). Then 
Pa- satisfies the axial symmetric wave equation 

d^tP„{z,r,t) = {r-^drrdr + dl)P,{z,r,t), (z, r, t) G K x (0, cx))^ , (2.7) 

P,(z,r,0) = F,(z,r), (z, r) G M x (0, c») , (2.8) 

9tP<,(z,r,0) = 0, (z,r) e M X (0,cx)) . (2.9) 

Moreover P„ remains bounded as r ^ 0. 
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Proof. Equations (2.8), (2.9) and the boundedness as r ^ immediately follow from (2.2), (2.3) 
and the definitions of and F^, . In the cylindrical coordinates , the Laplace operator is given 
by the well known expression 

A = r-^drrdr + d^, + r-^dl. 
Integrating the equation Ap ~ d^p with respect to a yields (2.7). □ 
Note that (2.7)-(2.9) is uniquely solvable if we require that its solution remains bounded as 

Remark 2.2. Proposition 2.1 is the basis of the following two-stage procedure which reconstructs 
the initial pressure f in (2.1)-(2.3) from the data {Ga)a-eS^' 

(i) For a G (fixed position of the stack of circles) determine the initial pressure Fa- of 
(2.7)-(2.9) using data G^. 

Repeating this procedure for every a, one obtains a family of functions F^, a G , corre- 
sponding to averages over circles centered on d(Bii(0) x R). 

(a) Next one recognizes that for fixed z — zq, the function 

{a,r) Fa{zo,r) 

is the circular mean transform of f\^z=za} with centers on a circle. For the circular mean 
transform stable analytic inversion formulas have been discovered recently [2, 8, 11, 18, 19]. 
Exemplarily, one of the formulas in [8, Theorem 1.1] reads 

/(x', zo) = ^ (^j^\drrdrFa){za,r) log ^ - |x' - Ra\'\ drj da , (2.10) 

where x' denotes the coordinates in the plane x {zq}. 

The key task for reconstruction / is to derive stable and fast algorithms to reconstruct the 
initial data in (2.7)-(2.9) from measurement data Go-. A possible reconstruction method is based 
of time reversal (back-propagation) similar to [4, 5, 14]. However, the degeneration of r~^drrdr at 
r = and the open detector set may cause difficulties in such procedures. The inversion approach 
in this paper is based on analytic inversion formulas for reconstructing F^. 

2.2 Exact Inversion Formula 

In the following denote by 



S{0}(w):=/ (j){t)sin{ujt)dt, 0e ii((0,oo)),w > 0, 
Jo 

/•oo 

C{0}(w):=/ <f){t) coa{ujt)dt , e i^((0, oo)) , w > , 



F{0}(fc):= / 0(z)e-''''^dz , (/) e L\R) , k e R , 
Jr 

/•oo 

n{(f>}{v):= <f){r)Jo{vr) rdr , (/) € L\{0,oo),r^/^dr) ,v > , 
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the sine, cosine, Fourier, and Hankcl transform, respectively. (Here Jq is the zero order Bessel 
function [1].) When the above transforms are apphed to functions depending on several variables 
then the transformed variable is added as subscript, e.g., H,. {i^o-} {z, v) = F„{z^ r)Jo{vr)rdr. 



Proposition 2.3. Let f e Cf^ {Br{0) xR) and define F^{z,r) andG„{z,t) hy (2.5), (2.6). Then 
the relation 



H..(F.i^.i)(M) = ^^''";/°'"';;„-^' (2.11) 

holds whenever Jo(?'ciotiO 7^ 0. 

Proof. The zero order Bessel function satisfies r~^drrdrJo{r) = —Jo{r) on (0,oo). Hence the 
chain rule implies that r~^drrdrJQ{rv) ~ — u^Jo(ru) for every r,v > 0. Separation of variables 
shows that the functions 



{z,r,t) 1-^ e*''^ cos(tVfcM^) Jo(7-y), {k,v) € R x (0,oo) , 

solve (2.7), (2.9). Employing the initial condition (2.8) and the inversion formulas for the Fourier 
and Hankcl transforms it follows that the unique bounded solution of (2.7)-(2.9) is given by 

P„{r,z,t) := ^ [ I F{k, v)e''''' Jo{rv) cos{t^ + v^)vdvdk , (2.12) 
Jk jo 

with F{k,v) := U,{F , {F„}} {k,v). 

Substituting w = -y/fc^ + in (2.12) and putting r =■ rdct afterwards, leads to 

G„{z,t) = ^ J (^J^ JoirdetVuJ^ - k^)LoF{k, y^uj^ - k^) cos(wt)dwj e'''^dk (2.13) 

The inversion formulas for the Fourier and Cosine transforms now imply that 

r< f-i? w fi, f Jo(rdot\/tj2~~P')cjF(/c, \/w2~~fc2), if > fc , 

Q{F.{G4}(fc,a;)--| otherwise. 

Solving the last equation for F shows (2.11). □ 
Proposition 2.3 implies that F^ can be reconstructed from data Go- as follows: 

(i) The data are Fourier and cosine transformed, yielding to {F^ {Go-}}. 

(ii) According to (2.11), Ct {F, {G^}} is mapped to {F;, {F^}}. 

(iii) Finally, application of the inverse Fourier and Hankcl transforms yields 

F,{z, r) = -!- / / {F, {F,}} (fc, v)Jo{rv)e'''' vdvdk . 
Jk Jo 

Remark 2.4 (Instability of (2.11)). Inversion formula (2.11) is not defined when Jo(f'dotw) equals 
0. From the proof of the above theorem it is clear that for exact data 

Ct{F4Go}}(/c,y^;iTF) =0, neN, (2.14) 
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Figure 2: Cross section of experimental buildup when rdot < R (left) and rdct > 2i? (right). In 
both case object has to be supported in the gray disc. 



with (w,i),igN denoting the zeros of r i~> Jo{r^ctv). In practice, however, only noisy (approximately 
measured) data ~ Ga are available. In general, 

Q{F.{G^}}(fc,y^;jTI^) ^0. 

It is therefore difficult to stably evaluate the quotient in (2.11) in practice. 

Because the acoustic pressure is measured outside of the investigated object, only the following 
two situations occur in practical applications, see Figure 2: 

(i) The stack of circles is strictly outside the object. In this case rdot < R and supp(/) C 
B,j„,,,,(0) xR. 

(ii) The object is enclosed in the stack of circles. In this case rdct > 2i?. 

For the case rdct > 2i? we will provide two stable formulas based on expansions in bases of special 
functions. The case rdct < R turns out to be harder, we currently do not have a stable alterative 
to (2.11). In this case the function F {Fo-} (fc, •) is not supported in the interval (0, rdct) and thus 
it cannot be expanded into a Fourier Besscl series which is crucial in the proves of theorems 3.1 
and 3.4. However, in the limiting case rdct <C i? a stable reconstruction formula is obtained, 
compare with Remark 3.6 below. 

3 Stable Inversion Formulas 

In the following we fix / S C^{Br{Q) x R), ct e S"!, define F„, G„ by (2.5), (2.6), and let (w„)„eN 
denote the zeros of the function v i— > Jo (rdct "y)- 

Our first stable inversion formula is as follows: 
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Theorem 3.1 (The D'Hospital trick). Assume r^ct > 2i?. Then 
2 



' det 



E 



St{tF,{G„}}{k,^/k^T^)v^ Jo{rv^) 



e^'^'^dk 



(3.1) 



for any (z,r) e R x (0,oo). 



Proof. The assumptions / € C^(i?fl(0) x M) and rdct > 2_R imply that Fz{Fcr}(fc, •) is compactly 
supported in (0, rdot)- It can therefore be expanded in a Fourier Bessel series [26] 



{F„} (k, r)^4-Y. {F. {F,}} (k, v„) /'^''''"\, 



' det 



According to (2.11) we have 



iIr{F,{F,}} {k,v) 



2Ct{¥,{G,}}{k,^/WTv^ 



Applying the rule of D'Hospital gives 



V ^ { i;„ : n e N} 



. .... ^ 2 d/dv \Ct {F, {GA} (k, Vk^ + v^)] 
llr{F,{F„}}{k,Vn) = - lim — ^ \ ^ \ ' ^ 

TT v^v„ g/Qy [Joirdotv)V k'^ + D^J 

2 St{tFAG^}}{k,V^T^)^^ 
= — lim - — - — 



^ 2 S,{tF,{G^}}(fc, VP+t;2)«„ 

Trrdet Jl{rdetVn){k'^ + Vn) 

Inserting (3.3) in (3.2) and using the Fourier inversion formula shows (3.1). 
Remark 3.2. [Stabilty of (3.1)] From the asymptotic approximation (see [1]) 

2 / TOTT TT 



Jra ('^) 



2 4 



^ , for X oo , 



of the m-th order Bessel function it follows that 

Vn^ , \Jl{rdctVn)\ - 

rdct V TrrdotWn 

Moreover the summands in (3.1) take the asymptotic form 



for n — > oo . 



St {tF,{G^}} (k^^^T^) 



k^ + Jl{rdetVnY 

1 |cos (i;„-f)|(2/(7rr^„))V^ 



< 



St{tF,{G.}}{k,./WT^, 
'''^^ St{tF,{G^}} {k,^k^+vl) 



4r 



3/2 
det 



(3.2) 



(3.3) 
□ 



Consequently, the parts that do not depend on the data G^ are bounded, and (3.1) can be imple- 
mented in stable way. 
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In the sequel we derive an additional inversion formula that circumvents the division by zero 
problem. In fact, our formula will be a consequence of the following result derived in [30]. 

Proposition 3.3. Let p denote the unique solution of (2.1)-(2.3) and let /™ and denote the 
Fourier eoefficients of f{^a{z,r,a)) and 



g„{a,z,t) 

with respect to a. Then 



p{^aiz,r,a),t), fort>0, 
0, otherwise , 



H.(F.(/n}(M) = ^^^%|^ai%^S, (...).Mx(0.»), (3.4, 

with denoting the m-th order second kind Hankel function. 

Now the second stable inversion formula can be stated as follows: 
Theorem 3.4. Assume rjct > 2i?. Then 



Here Ga is extended by G^iz, t) ~ for t < 0. 

Proof. We use again the Fourier-Bessel series (3.2) of proof of Theorem 3.1. Recalling the defini- 
tions of Fcr, Ga and the Fourier coefficients /™, 17™ one notices that ~ f^, Ga = Therefore 
(3.4) for TO = implies 

H.{F.{F4}(fc,.) ^ 2FaF {G4}(fe,v/W)^ ^^^^^ ^ MX (0,00). (3.6) 

Inserting (3.6) in (3.2) and using the Fourier inversion formula shows (3.5). □ 

Equation (3.6) is quite similar to (2.11). However, in the denominator in (3.6) the zero order 
second kind Hankel function appears (instead of the the zero order Bessel function) which cannot 
be zero for a finite argument [1]. Moreover, the asymptotic expansion of the Bessel and the 
second kind Hankel function show that the summands in (3.5) that do not depend on the data 
Ga remain bounded as n — > 00. 

Remark 3.5. The derivation of (3.4) is based on the following Green function expansion in 
cylindrical coordinates [30] 

-iaj|*„(z,r,a)-*„(zo,rdct,ao)l 



\^aiz,r, a) - ^a{zo,rdct,ao)\ 

■JIYI ^™(^'^' ^^dot)e-™("-"'') J e-'-(--o) dz , (3.7) 



'ITT 
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with V = sign.{uJ)^/\uJ^^^W\, 



{vrdet)J7n{vr), 



if uj^ > k"^ , 



2i/7TKmi\v\ rdet)l7ni\v\ r), otherwise, 



and Irn, Km denoting the m-th order modified Bessel functions of first and second kind, respec- 
tively. Here its is assumed that r^ot > 

Interchanging the roles of r and r^^t implies that for r^ot > f the Green function expansion 
(3.7) holds with 



Ara{vr,vrdct) = \ r 7| I \J /: | > 



therwise . 

Similar to the proof of (3.4) in [30] this leads to a formula for reconstructing F„ in the case 
Tdct < 2i?, however again an unstable one with Jo(i'fdct) in the denominator. 




5 10 15 20 25 30 35 40 




0.2 0.4 0.6 O.i 



1.2 1.4 1.6 



Figure 3: The first 50 denominators MrdctVn/K) in (3.10) for K = -i (left) and K = 100 (right). 



Remark 3.6. Suppose that r^ct < R and that f is supported in Bn-r^^^ (0) x M. For ri > 2R—rdct 
let (un)neN denote the zeros of v t-^ jQ(riv). Then one can expand Fz {F^} (k, ■) in a Fourier 
Bessel series (see [26]) 

Fz{F^}{k,r) = ^Y.^r{FAFa}}{k,Vn)4T^^ (fc, r) G M x (0, oo) . (3.8) 

According to (2.11) we have 



H,(F.(f,H(..., = i5i£ii^Mi(|^3 „,K:„eNl. (3.9, 

Jo(fdotf)vfc-^ + 

If we assume thatri is a integer multiple ofr^ct, i-e., ri = Kr^ct, thenvn = Vn/K ^ {v„i : m e N} 
for any n eN. Therefore, inserting (3.9) in (3.8) yields 
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In general, (3.10) is a again sensible to noise when Vn gets close to a zero o/ Jo(rdotw). 

In the limiting case rdot ^ R o-nd for n not to large, the denominators Jo(riu„) are well 
bounded from below (see right image in Figure 3). In this case, truncating (3.10) leads to a stable 
inversion formula. In particular, for rjot = one obtains 



/or any ri > 2R. Together with (2.10) this provides a novel inversion formula for PAT using 
point-like detectors on a cylindrical recording surface. 

4 Numerical Experiments 

In practice one deals with discrete measurement data 

Gi[m,ii] := G^i(z„,<n), (l,m,n) £ {1,...,N^} x {l,...,Nj x {l,...,Nt}, 

where G„ is as in (2.6), and where CTi = 27r(l — l)/No., = H{m — l)/^^ and = T(n — 1)/Nt 
are discrete samples of the angle, height and time, respectively. Here H > represents the finite 
height of the stack of circular integrating detectors (see Figure 1) and T is such that Ga{z, t) = 
for t>T and z e [0,H]. 

In this section we outline how to implement (3.1) and (3.5) in order to find an approximation 

Fi[m,n] ~ F^,( ), (m, n) e {1, . . . , N^} X {1, . . . , Nr} , 

with 7'ii — r(jet(ii ^ 1)/Nr- Having calculated such an approximation, one can reconstruct a 
discrete approximation to / by applying the filtered back-projection algorithm of [8] for fixed m, 
see Remark 2.2. 

A numerical reconstruction method based on (3.1) is as follows: 

(i) The discrete Fourier transform (with respect to the first component) of the data 

F{Gi} [m,n] := ^ GiK,n] e-*2™(m'-i)/N, ^4_;^) 

m' = l 

with (m,n) G {— Nz/2, . . . , N2/2 — 1} x {l,...,Nt}, is considered as an approximation to 
F{F,J(27rm/ff,t,). 

(ii) The sine transform S {i F {i^o-i}}, evaluated at 



Wm,n = V(27nn/^)' + ^'n,(m,n) G {-N^/2, . . . , N,/2 - 1} x {0, . . . , - 1} , 
is approximated by the trapezoidal rule, leading to 

Nt 

S {tF {Gi}} [m, n] ^ 4'F {Gj [m, n'] sin(w„,nin') ■ (4.2) 

n' = l 

(iii) Finally, truncating the Fourier Bessel Series and approximating the inverse Fourier trans- 
form with the trapezoidal rule leads to discrete version of (3.1): 

V Irr, r.]- '^^ "'-^ ^ "^^^ «n'S{tF{Gi}} [m',n'] ^ -«27rm'(m-l)/N, 

with (m, n) e {-N^/2, . . . , Nj2 - 1} x {0, . . . , N,. - 1} in formula (4.3) 
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A numerical reconstruction method using (3.1) is be obtained in an analogous manner. In 
this case one replaces (4.3) by 

T7 r,l • 4T "'.^ ^ F {F {Gi}} [m',n'] n -j2™'(m-i)/N, fA 

Fim,n.= — X F7T7 TTT Jo(?'nWn')e ^ " , (4.4) 

'^MotNt ^i-'o W„,/,n'i?o (''dctt'n')'^o(''dctVn')^ 

which is the discrete analogue of (3.5). 

To give a rough estimate of the computational complexity for the previous calculations let us 
assume = = Nt = =: N and that the values of the sine function and the Bessel function 
are pre-computed and stored in lookup tables. Then the evaluation (4.1) needs ©(N^ logN) float- 
ing point operations (FLOPS) whereas (4.2) and (4.3) require 0(n3) FLOPS. The filtered back 
projection formula (2.10) also requires 0{'N^) FLOPS, see [8]. For three dimensional reconstruc- 
tion (4.1), (4.2), (4.3) and the filtered back-projection formula have to be applied N times. Hence 
the total number of FLOPS is estimated as 

Nflops = N (0(n2 logN) + 0(n3) + 0(N^)) = 0(n4) . (4.5) 

Note that three dimensional back-projection type formulas which use point measurement data 
have complexity ©(N^). 

In the following numerical experiments we take R = 0.4, H = 3.75 and T — A. The synthetic 
initial data / is assumed to be a superposition of radially symmetric objects around centers x„, 
i.e., 

/(x) = ^/„(||x-x„||), xeM^^ 

n 

The acoustic pressure generated by a single radially symmetric object at position x and time t is 
given by (sec [13]) 

Pn(x, t) = "^^^'^'l^'^^ /„ (I ||x - x„|| - t\) . (4.6) 
By the superposition principle the total pressure is 

N 

p{x,t)^"^p^{x,t), (x,t) e X (0,c^). 

n=l 

The measurement data Gcr{z,t) = l/(27r) J^^ p{^aiz,rdot,ce),t)da, see (2.4), (2.6), were gen- 
erated by evaluating of (4.6) followed by numerical integration over a. The radius rdot of the 
circular integrating detectors is chosen to be 2R. In this case the stack of circular integrating 
detectors fully encloses the synthetic initial data /, see right image in in Figure 2. 

Figure 4 shows a vertical cross section of the initial pressure / and the measurement data Go- 
where Gaussian noise with a variance of 10% of the maximal data valued is added. The stack 
of circular integrating detectors is centered to the left of the objects. The presented discrete 
implementation Nt = 320 measurements in time and = 300 in space. In both reconstructions 
the value was chosen to be = 130 

The reconstructions of Fo- with (4.3) from exact and noisy data are depicted in Figure 5 from 
formula (4.3) and with formula (4.4) in Figure 6. In the reconstructed images one notices some 
blurred boundaries which are limited data artifacts [20, 24, 29] arising from the finite height of 
the stack of circular integrating detectors. Moreover the images reconstructed with (4.4) are less 
sensitive to noise. 
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Figure 4: Left: Cross section of five absorbing spheres {z versus r). Right: The measurement 
data with 10% Gaussian noise added (z versus t). 



5 Conclusion 

In this article a novel experimental buildup for PAT using circular integrating detectors was 
proposed. For collecting measurement data a fiber-based Mach-Zehnder or Fabry-Perot interfer- 
ometer can be used as an circular integrating detector. We showed that the 3D imaging problem 
reduces to a series of 2D problems. This decomposition can be used to reduce the operation count 
of derived reconstruction algorithms. We derived two stable exact reconstruction formulas, (3.1) 
and (3.5), for the case that the object is contained in the stack of detecting circles. In the case 
where the object is outside the detecting circles, a stable reconstruction formula is obtained for 
the limiting case rdct R- As a byproduct, this leads to a novel reconstruction formula (3.11) 
for PAT using point detectors on a cylindrical surface. 
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Figure 5: Reconstruction with (3.1) from simulated (left) and noisy data (right). 
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